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Abstract 

The Kohn-Sham equation is a powerful, widely used approach for computation of ground 
state electronic energies and densities in chemistry, materials science, biology, and nanosciences. 
In this paper, we study the adaptive finite element approximations for the Kohn-Sham model. 
Based on the residual type a posteriori error estimators proposed in this paper, we introduce 
an adaptive finite element algorithm with a quite general marking strategy and prove the 
convergence of the adaptive finite element approximations. Using Dorfler's marking strategy, 
we then get the convergence rate and quasi-optimal complexity. We also carry out several 
typical numerical experiments that not only support our theory, but also show the robustness 
and efficiency of the adaptive finite element computations in electronic structure calculations. 
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1 Introduction 

The Kohn-Sham density functional theory (DFT) [35, 36. is a powerful, widely used approach for 
computation of ground state electronic energies and densities in chemistry, materials science, biol- 
ogy, and nanosciences. Consider a molecular system consisting of M nuclei of charges {Z\, • • • , Zm} 
located at the positions {Ri, • • • , R-a/} and N electrons in the non-relativistic and spin-unpolarized 
setting. The ground state solutions of the system can be obtained by minimizing the Kohn-Sham 
energy functional 
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with respect to the Kohn-Sham orbitals {<fo}£Li under the orthogonality constraints 

j>i<t>j = Sij, 1 < i, j < N, 



M 



z k 



where V e xt(%) = ~ 1 "73 — r ^ s the electrostatic potential generated by the nuclei, p(x) 

JV 

|0i(a;)| 2 is the electron density, and e xc (p) denotes the exchange-correlation energy per unit 

i=l 

volume in an electron gas with density p. 

Since the core electrons do not participate in the chemical binding and remain almost unchanged, 
a pseudopotential approximation is usually resorted to in practical computations of the Kohn-Sham 
model, which is to replace the strong Coulomb potential of the nucleus and the effects of the core 
electrons by an effective ionic potential acting on the valence electrons. Therefore, under the 
pseudopotential framework, only valence electrons are involved. The pseudopotential consists of 
two terms: a local component V\ oc (whose associated operator is the multiplication by the function 
Vioc) and a nonlocal component V n \ (an operator whose expression is given in Section [2]). As a 
result, the energy functional is replaced by 



N N 

E(mli) = TjE / \VMx)\ 2 dx+ / V loc (x)p(x)dx +J2 
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where N now denotes the number of valence electrons and {4>i}iLi denotes the set of the pseudo- 
orbitals of the valence electrons. Note that if we take V\ oc — V cx t and V n \ = 0, then (| 1 . 2() becomes 
(jl.ip . Consequently, we restrict our discussions to minimizing energy functional (|1.2[) in this paper. 

The Euler-Lagrange equation corresponding to this minimization problem is the so-called Kohn- 
Sham equation: 

{-\A + V eS {p))4>i = Hi<t>i inM 3 , i= 1,2, ••• ,N, 

(1.3) 

bj = Sij, 



where /Xj (i — 1,2, •■ • , N) are N lowest eigenvalues, V s(p) is the effective potential relative to 
the last four terms in energy functional (11.21) . This is a nonlinear integro-diffcrcntial eigenvalue 
problem, which is also called self-consistent field (SCF) equation as to emphasize the nonlinear 
feature encoded in V e s(p)- 

We understand that the Kohn-Sham approach achieves so far the best balance between accuracy 
and efficiency among all the different formalisms of electronic structure theory, and simulations of 
large-scale material systems with Kohn-Sham DFT are still computationally very demanding (say, 
thousands of electrons or more). As a result, efficient numerical algorithms that can be scalable on 
parallel computing platforms are desirable to enable DFT calculations at larger scale and for more 
complex systems. We see that real-space techniques and methods for electronic structure calcula- 
tions have been derived much attention from scientific and engineering computing communities and 
remarkably developed during the last two decades, among which the finite element method possesses 
several significant advantages [5J [251 ESI EH [SHI [S7] . Although the finite element method uses more 
degrees of freedom than that of traditional methods like plane waves and Gaussians, strictly local 
basis functions produce well structured sparse Hamiltonian matrices; arbitrary boundary conditions 
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can be easily incorporated; it is scalable on parallel computing platforms; more importantly, it is 
variational and relatively straightforward to implement adaptive refinement techniques for describ- 
ing regions around nuclei or chemical bonds where the electron density varies rapidly, while treating 
the other zones with a coarser description, by which computational accuracy and efficiency can be 
well controlled. 

We observe that even in the pseudopotential setting, the eigenfunctions of (|f .31) still vary rapidly 
around nuclei or chemical bonds [BJ 1181 132] . Hence it is very natural to apply adaptive finite ele- 
ment (AFE) approaches to improve the approximation accuracy and reduce the computational cost. 
Indeed, we see that AFE computations have been quite successfully used in solving Kohn-Sham 
equations and electronic structure calculations. Tsuchida and Tsukada combined the finite ele- 
ment method with the adaptive curvilinear coordinate approach for electronic structure calculation 
of some molecules [5FJ [59] ; Shen and Zhang introduced some adaptive tetrahedral finite element 
dicretizations in their theses [511 W2\ and calculated several typical molecular systems efficiently 
[3"21 1521 [Ml R)3"] ; Bylaska et.al used adaptive piecewise linear finite element method on completely 
unstructured simplex meshes to resolve the the rapid variation electronic wave functions around 
atomic nuclei [10 ; Dai et.al designed some parallel adaptive and localization based finite element 
algorithms for typical quantum chemistry and nanometer material computations containing more 
than one thousand atoms using tens of hundreds of processors on computer cluster (TTJ [T51 UHl HI] ; 
Gavini et.al constructed a finite element mesh using unstructured coarse-graining technique and 
computed materials systems [HI [55J ; Yang successfully scaled their AFE simulations to over 6000 
CPU cores on the Tianhe-IA supercomputer in his thesis [65] . The AFE simulations carried out in 
this paper also show the robustness and efficiency of the AFE computations in electronic structure 
calculations. We may refer to [571 |5S] and references cited therein for other interesting discussions 
on adaptive finite element method (AFEM). 

We see that it is significant to understand the mechanism of AFE computations, analyze the 
AFE approximations of Kohn-Sham equations, and give a mathematical justification of the AFE 
algorithm. We note that the AFE computations are based on some a posteriori error estimators 
and there are very limited work concerning analysis of the a posteriori error estimators and con- 
vergence of AFE approximations for DFT. In [TH [15], the authors of this paper considered the 
nonlinear eigenvalue problems derived from the orbital-free DFT and obtained the convergence and 
optimal complexity of the AFE algorithm. We understand that the orbital-free DFT is viewed as 
a simplification of the Kohn-Sham DFT, in which only one eigenpair is involved. In this paper, 
we shall propose and analyze two AFE algorithms for Kohn-Sham DFT calculations and study the 
associated convergence and optimal complexity. Indeed, this is one of a series of work of our group 
on AFE computations for electronic structures. 

Let us now give a somewhat informal description of the main results of this paper. We propose 
and analyze two AFE algorithms: Algorithm [3J] and Algorithm |4T] which are based on the residual 
type a posteriori error estimators. We show the a posterior error estimates (see Theorem 14.3ft and 
prove that 

• Under some reasonable assumptions, all limit points of the AFE approximations of the ground 
state solutions are ground state solutions (see Theorem 13. 1|) . 

• Under other reasonable assumptions, some eigenpairs (in particular, ground state solutions) 
can be well approximated by AFE approximations with some convergence rate (see Theorem 
Ol) . 

We also study quasi-optimal complexity of AFE approximations (see Theorem 14. 6[) . 

We mention that Algorithm 13.11 and Algorithm 14.11 may be viewed as some extensions of asso- 
ciated existing algorithms for linear elliptic partial differential equations of second order and have 
been in fact used for years, for instance, in package RealSPACES (Real Space Parallel Adaptive 
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Calculation of Electronic Structure) of the State Key Laboratory of Scientific and Engineering Com- 
puting, Chinese Academy of Sciences. As we see, the numerical analysis for AFE approximation 
is important and has been also derived much attention from the mathematical community. Since 
Babuska and Vogelius [3] gave an analysis of an AFEM for linear symmetric elliptic problems in 
one dimension, there has been much investigation on the convergence and complexity of AFEMs 
in literature (see, e.g., [T21 [23 1231 1301 [S3] and the references cited therein). In the context of 
the finite element approximations of linear eigenvalue problems, in particular, we see that there 
are a number of works concerning a posteriori error estimates [SJ HH [231 |331 133 GUI ED] , AFEM 
convergence [23 [23 1301 133 133] and complexity QS [23 [HI [33] • 

However, there are several crucial difficulties in numerical analysis of the Kohn-Sham equation: it 
is a nonlinear eigenvalue problem whose eigenvalues may be degenerate, and a number of eigenpairs 
must be involved; the associated energy functional is nonconvex with respect to density p 1 as a result, 
there is no uniqueness result for the ground state solutions; the energy functional is invariance under 
unitary transforms, which also induces redundancy of the ground state solutions. To handle these 
difficulties arise from the Kohn-Sham equations, we shall present some sophisticated arguments 
and consider the convergence under the distance between solution sets; investigate the convergence 
rate and optimal complexity under certain inf-sup assumption; and exploit the relationship between 
the finite element nonlinear eigenvalue approximations and the associated finite element boundary 
value approximations. 

The rest of this paper is organized as follows. In Section [21 we provide some preliminaries for 
Kohn-Sham DFT problem setting and residual type a posteriori error estimator based AFE meth- 
ods. We prove the convergence of AFE approximations in Section [3] and analyze the convergence 
rate and optimal complexity of an AFE algorithm in Section UJ In Section [5j we present some 
numerical experiments that support the theory. Finally, we give some concluding remarks. 

2 Preliminaries 

Physically, the Kohn-Sham equation is set in R 3 . However, due to the exponential decay of the 
ground state wavefunction of the Schrddinger equation (c.f., e.g., [21 [ST]) and the fact that Kohn- 
Sham model is an approximation of Schrddinger equation, R 3 is usually replaced by some polyhedral 
bounded domain Q, C R 3 in practical computations for Kohn-Sham equation. 

For k £ M. NxN , we denote its Frobenius norm by \k\. For p > 1 and s > 0, we denote by W S ' P {VL) 
the standard Sobolev spaces with the induced norm || • || s ,p,n (see, e.g. [3 US]). For p = 2, we denote 



by H s (n) = W S < 2 {VL) with the norm || • || s , n = || • || s , 2 ,n, and = {v e H 1 ^) : v \ dn = 0}, 



where v |an= is understood in the sense of trace. The space H 1 (fi), the dual of i?g(fi), will also 
be used. Let T~L = (Hq(£1)) n be the Hilbert space with inner product 



($,*) = ^/ 0,-0, for $ = (</>!,••■ ,<?W)>* = Oi,-" ,iPn) G H. 




Let Q be a subspace with orthonormality constraints: 



Q = {$ G H : = I NxN }, 




For $ € H and a subdomain u C U, we shall denote by 
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N 

p$ = |0j| 2 an d (sometimes abuse the notations for simplicity) by 

i=l 

/ N \ 1 / 2 / N \ Vf 

ll^lko, = f S J = 0,1; ||$llo,p )W = (Ell^HS*, J ; !<P<6- (2.1) 

In our discussions, we shall use the following sets: 

S NXN = {M g R JVxAr . M T = M}) = {M g R JVxAT . # T = _ M} ^ 

We may decompose H into a direct sum of three subspaces (see, e.g., [IS]): 

n = 5$ © A* © 7i 

for any $ G Q, where 5$ = <S>S NxN , A<s, = <S>A NxN , and T$ = {f£«: = G R NxN } . 

For convenience, the symbol < will be used throughout this paper, and A < B means that 
A < CB for some constant C that is independent of mesh parameters. We use ^(p, {c\,Ci)) to 
denote a class of functions satisfying some growth conditions: 

&>{p, (ci,c 2 )) = {/: 3 a X) a 2 G R such that c x t v + at < f{t) < c 2 t p + a 2 V t > 0} 

with ci G K and C2,p G [0, oo). 

2.1 Problem setting 

We consider the following general form of Kohn-Sham energy functional 



r (i N N \ l 

^ n V i=l i=l / 



P*) (2-2) 



for $ = (0i, ^2, • • • , 4>n) G which includes the cases of Coulomb potentials and pseudopotential 
approximations. For the Coulomb potential setting, V\ oc = — J2k=i \x-R k \ an< ^ = 0- While 
for pseudopotential approximations, V\ oc is the local part of pseudopotential and Vni is a nonlocal 
pseudopotential operator (see, e.g., [40]) given by 

n 

K!0 = ]>>,G)G 

3=1 

with Q G L 2 (Q)(j = 1, 2, • • • , n). D(p$, p$) is the electron-electron Coulomb energy defined by 

D(f,g)= [ f(g*r- 1 )= [ [ f( x )g( y )-^—dxdy, 
Jn JnJn \ x ~ V\ 

and e xc (t) is some real function over [0, oo). In our following analysis, we require V\ oc belongs 
to L 2 (tt). We point out that V\ oc G L 2 (ty is a very weak condition, which is satisfied by both 
the Coulomb potential Vext(x) = — Ylk=i la-RJ an< ^ ^ c l° ca l P ai 't of pseudopotential. Since 
e xc : [0, oo) — > R does not have a simple analytical expression, we shall use some approximations 
and assume throughout this paper that 

e xc (t) G 5»(3, (ci, c 2 )) with ci > or e xc (i) G £»(4/3, (ci, c 2 )), (2.3) 
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which is satisfied by almost all the LDAs. 

The ground state of the system is obtained by solving the minimization problem 

inf {£($):$ G Q}, (2.4) 

and we refer to [3J [TTJ [TB] for the existence of a minimizer. Note that the energy functional (|2.2p is 
invariant with respect to any unitary transform, i.e. 

N 

= E(9U) = £(£^0 V £7 = K)5 =1 G M , (2.5) 
j'=i 

where Q NxN is the set of orthogonal matrices. It follows from (j2.5[) that if $ is a minimizer of 
l|2.4p . then $£/ is also a minimizer for any orthogonal matrix U. For any $ £ ?{, we define the 
equivalence class 

[*] = {9U, V U G G A,xAr }. 

We see that any minimizer $ = • • • , 0jv) of (|2.4|) satisfies the following weak form (i.e. the 
Euler-Lagrange equation associated with the minimization problem): 

JV 

(fr»0j,«) = (^A^u) VuGi^Q), i= 1,2,- ■■ ,A, 

j=i (2.6) 



where H$ is the Kohn-Sham Hamiltonian operator as 



= ~A + V Xoc + V ai + [ j^\dy + e'^M (2.7) 
2 Jn\- -y\ 



and 



A = M?j=i = ( I 

vn / 1.1=1 



is the Lagrange multiplier. Since the uniqueness of the ground state solution is unknown even up 
to a unitary transform, we define the set of ground states by 

9 = j(A,$) G R NxN x Q : E(<$>) = mmE(^) and (A,$) solves (|2T5])| . (2.9) 

Note that the electron density /?$ and the operator H$> are also invariant under any unitary trans- 
form, we may diagonalize the matrix of Lagrange multipliers A. More precisely, there exists a 
U G NxN , such that the Lagrange multiplier is diagonal for ^ = <&U = (ipi, ■ ■ ■ i.e., 

Consequently, instead of (|2.6[) , we may consider a form with diagonal multiplier as follows: 

= VweflJcn), * = 1,2,- •• ,iv, 

/ / A (2 - 10) 



G 



which is the well known Kohn-Sham equation. 

Note that any solution of (12.61) can be obtained from a unitary transform of some solution of 
(|2.10p . That is, if * is solution of (|2.10p . then any $ G [*] is a solution of (|2.6|) . Consequently, we 
also call (|2.6[) Kohn-Sham equation. 

It is well known that the ground state has one electron in each of the N orbitals with the lowest 
N eigenvalues. Therefore, the ground state solutions in (|2.9[) can be obtained by solving the lowest 
N eigenpairs of (|2.10[) . 

For convenience, define T : R NxN x H — >• H* by 

N N 

(j-(A, $), r) = J2 (H*<Pi - E x a<t>hn) v r = ( 7i )£i e n. 

i=l 3 =1 

The Frechet derivative of T with respect to $ at (A, <!>) is denoted by J"4(A, $) : H -> H*, where 

1 W 

(j-4(A,$)*,r) = ^"($)(*,r) - £ (A«^,7i) 

N N N N 

i=l j=l ij— 1 i,j'=l 

To study the convergence and complexity, we need the following assumptions 
Al |e' c (i)| + |t«&(i)| G ^(pi, (ci,c 2 )) for somepi G [0,2]. 

A2 There exists a constant a G (0, 1] such that |e" c (£)| + |*e^(t)| < 1 + t Q_1 V i > 0. 

A3 (A, $) is a solution of (|2.6[) and there exists a constant /? > depending on (A, $) such that 

i.l,»2>* (2.11) 
rer** e7i ||* x^i r||i, n 



Remark 2.1. Assumption Al ensures some lower semicontinuity property of the energy functional. 
Assumption A2 is used to make sure the continuity of J-"4(A, $). We see t/tat Assumption A2 
implies Assumption Al and i/ie commonly used X a and LDA exchange-correlation energy functional 
satisfy Assumption A2. 

Assumption A3 is equivalent to that ^^(A, <£>) is an isomorphism from 7$ £o 7$. We observe 
that if Assumption A3 is satisfied for $ € Q, i/ien Assumption A3 is satisfied for any $ G [<&] wii/i 
i/ie same constant j3, too. We see that a stronger condition than (|2. 1 lj) i/iai 

(Ji(A,*)r,r) >7l|r||?, n vrer* 

is used in j7.?t 150) /. which is satisfied for a linear self-adjoint operator when there is a gap between 
the lowest Nth eigenvalue and (N + l)th eigenvalue \5(tf . 



2.2 Adaptive finite element approximations 

Let d n be the diameter of £1 and {Th} be a shape regular family of nested conforming meshes over 
fl with size h G (0, d n ): there exists a constant 7* such that 

^< 7 * VreTft, (2.12) 

Pt 
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where h T is the diameter of r for each r € 7/j, Pt is the diameter of the biggest ball contained in r, 
and ft, = ma.x{h T : t 6 7^}. Let denote the set of interior faces (edges or sides) of 7~h- 
Let S h ' k (Ti) be a subspace of continuous functions on CI such that 

S h - k (fl) = {ve C(fi) : »| T e P r fe V r e T4, 

where is the space of polynomials of degree no greater than k. Let S ' (CI) — S h ' k (Q) H Hq(CI). 
We shall denote Sg'^fi) by Sq(SY) for simplification of notation afterwards and let Vh = (Sq(CI)) n . 
We now consider the following finite element approximations of (|2.4[) : 

■w£{E(9> h ):$ h eV h n®}. (2.13) 

We refer to [51 [T3] for the existence of a minimizer of (|2.13l) . Note that any minimizer $/> = 
[4>i,h: ■ • • j 4>N,h) 01 (|2.13|) solves the Euler-Lagrange equation 

(H$ h <f>i, h ,v) = C^\ij,h<t>j,h,v) VuG fift(n), i = 1,2, ■■ ■ ,iV, 

i=i (2.14) 

k i,h<f>j,h = % 
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with the Lagrange multiplier 

A/, = (A lJ ,/ l )^ )=1 = / 4>i,hHq> h 4>ij L ) 
Define the set of finite dimensional ground state solutions: 
e h = UA h ,$ h )€R NxN x(®nV h ) :E($ h )= min and (A fc) solves (£H) ] ■ 

We have from [13] that the finite dimensional approximations are uniformly bounded, i.e., there 
exists a constant C such that 

sup (||4h||i in + |Afc|) <C. (2.15) 

{A h ,<s> h )ee h jie(o,d n ) 

Using a unitary transform, we can diagonalize and obtain a discrete Kohn-Sham equation 

(Hv h i/> iih ,v) = (m, h i/Ji,h,v) Vi)e5f(fl), i = 1,2,--- ,N, 

(2.16) 

^i,htpj,h = Sij 

with fii ih = (iff,,^,/,,^,/,). 

Similar to the continuous case, we also easily have that if "J^ is a solution of (|2.16|) . then any 
&h G [^h] is solution of (|2 . 14[) . Therefore, to get O/j, we always resort to solving (|2 . 16[) in practice. 

An adaptive mesh-refining algorithm usually consists of the following loop: 

Solve -» Estimate -> Mark -> Refine. 

Solve. This step computes a piecewise polynomial finite element approximation with respect to 
a given mesh. To simplify the analysis and do as the most work on numerical study of convergence 
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of AFE approximations, we shall assume throughout this paper that we have the exact solutions of 
discretized problems- 
Estimate. Given a partition Tt and the corresponding output (A^, from the "Solve" step, 
"Estimate" computes the a posteriori error estimator {i]h{$h) T )}r£T h i which is defined as follows. 
First, we define the element residual lZ T ($h) and the jump J e ($h) by 

N 

H T {<&h) = (H&^h - E ^ij,h<t>j,h) i=1 in r e %, 

3 = 1 

(\ 1 \ N 

where e is the common face of elements t\ and T2 with unit outward normals n\ and «2 , respectively. 
Let u>h{e) be the union of elements that share the face e, and uih(r) be the union of elements that 
share an edge with t. For r G Th, we define local error indicator 7y^($/j,r) and the oscillation 
osch($h,T) by 

r l l{<5> h ,T)=hl\\'R T {<S> h )\\l T + h e\\Je(®h)\\l, e , 



osc h (® h ,r) = h T \\n T ^ h ) - 7e T ($ ft )||o, T1 

where w is the i 2 -projection of w £ L 2 (Q) to polynomials of some degree on r or e. Given a subset 
uj a Q, we define the error estimator rjh w) and the oscillation osc/ l ( < l ) ? l , w) by 

r£ = ^ t) and osc£(* & ,w) = ^ t). 

Mark. We shall replace the subscript ft, (or /ij.) by an iteration counter k whenever convenient 
afterwards. Based on the a posteriori error indicators {r]k($k, T )}rGr fc > "Mark" gives a strategy to 
choose a subset of elements Aik of Tk for refinement. One of the most widely used marking strategy 
to enforce error reduction is the so-called Dorflcr strategy. 

Dorfler Strategy. Given a parameter < 9 < 1 : 

1. Construct a minimal subset Mk of Tk by selecting some elements in Tk such that 

E l'( $ ^)>^^ $ ^)' ( 2 - 17 ) 



2. Mark all the elements in Aik- 

A weaker strategy, which will be addressed as "Maximum Strategy" , only requires that the set of 
marked elements Mk contains at least one element of 71- holding the largest value estimator [2"9"II5U] . 
Namely, there exists at least one element T™ ax G Aik such that 

77 fe ($ fe ,T fc max ) = max77 fc ($ fc ,r). (2.18) 



1 Similar conclusion can be expected for the case where the errors of numerical integrations and nonlinear algebraic 
solvers are included (see Section [6)1. And we understand that the assumption is indeed a very important practical 
issue. 
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It is easy to check that the most commonly used marking strategies, e.g., Dorner's strategy and 
Equidistribution strategy, fulfill this condition. 

Refine. Given the partition Tk and the set of marked elements A4fc, "Refine" produces a new 
partition Tk+i by refining all elements in M. k at least one time. We restrict ourself to a shape- regular 
bisection for the refinement. Define 

K n ^T h+ i=T k \(T k nT k+ r) 

as the set of refined elements, we have A4 k c lZj- k ^j- k+1 . Note that more than the marked elements 
in Mk are refined in order to keep the mesh conforming. 



3 Convergence of adaptive finite element approximations 

In this section, we propose and investigate an AFE algorithm with Maximum Strategy for Kohn- 
Sham equations as follows: 

Algorithm 3.1. AFE algorithm with Maximum Strategy 

1. Pick an initial mesh To, and let k = 0. 

2. Solve (|2.16l) on Tk to get discrete solutions ipi,k){i = 1, • • • , N) and then Qk- 

3. Compute local error indictors rjk(^k, t) for all r G Tk- 
4- Construct Aik C Tk by Maximum Strategy. 

5. Refine Tk to get a new conforming mesh Tk+i ■ 

6. Let k = k + 1 and go to 2. 

We shall prove that all the limit points of the AFE approximations generated by Algorithm 13. II 
are ground state solutions of (|2.6[) . Given an initial mesh To, Algorithm 13.11 generates a sequence 
of meshes T\ , 75, ■ • • , and associated discrete subspaces 

S%°(a.) c s^(Q) c ... c s^(n) c s^ +1 (n) c ... c s^Q) c H^n), 

i #o (O) 

where Soo(fi) = U^S^fi) ■ Similar to the definition for V h , we set = (5' 00 (17)) JV . We 
have that Kx, is a Hilbcrt space with the inner product inherited from H and 

lim inf ||¥k-¥oo||i,n=(] V*ooel4o. (3.1) 

fe->°°* t e(s^(n))Jv 

Similar to (|2.4|) . we know the existence of a minimizer of energy functional (|2.2j) in n Q. 
We see that any minimizer $oo = (</>i,oo, ■ • • , 0at,oo) G Voo n Q solves the following Euler-Lagrange 
equation 

N 

{H 9ae ip i<00 ,v) = >Hj,ca<t>j,oo,v) V v e Soc{n), i = l,2,---,N, 

j=i (3.2) 

^.oo0j/,C3O — ^ij 



with the Lagrange multiplier 

Aoo = (Ay. oo)fj = i = ( / 4>j.ooH^> oo (f)i. 00 I . (3-3) 
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Define 

Oco - {(A^,^) € R NxN x (Voo HQ) : E^) = ^ min £(f) and (A^,^) solves 

Using similar arguments to those in the proof of Theorem 4.1 in [T3], we can prove that the AFE 
approximations converge to some limiting pair in 0oo. 

Lemma 3.1. Let {Ok}keN be the sequence obtained by Alaorithm \S.l\ We have 

lim E k = min E(V), 
lim d w (6 fc ,eoo) = 0, 

k— too 

where E k = £ , (<1>)((A, $) g Q k ) and the distance between sets X, Y C R NxN x H is defined by 

d n (X,Y)= sup inf (|A-jj| + ||$-*|| lin ). 

(A,$)ex UJ.*)ev 

Proof. Let (A fe , $ fc ) be the solution of (12TT4]) in R NxN x {V k C\Q) for fc = 1, 2, • • • , and {(A fcm , $ fe J} m£ N 
be any subsequence of {(Afc, <l>fc)}fc G N with 1 < k\ < k% < • • • < k m < ■ ■ ■ . 

Note that (|2.15[) and the Banach-Alaoglu Theorem yield that there exists a weakly convergent 
subsequence {<&k m . }jefi an d ^oo G Uoo H Q satisfying 

$fc mj ^ $oo in "H, (3.4) 

it is sufficient to prove 

£($00) = min E(V), (3.5) 

lim (||$ fcm . - + |A fem . — Aooj) = 0. (3.6) 

Since i?o(ri) is compactly imbedded in L p (il) for p e [2,6), by passing to a further subsequence, 
we have that $ km . — > $00 strongly in (L p (fl)) N as j — > 00. Thus we can derive from (|2.3[) that 

lim / ex C (p* fc ) = / e xc (p #( J, 

lim D(p^ k ,p^ h ) =-D(p*oc)P*oc) ! 

J—>00 3 3 



and hence 



liminf£($ fcm .) > £(*«,). (3.7) 



Note that p.ip implies that {<&fc m } is a minimizing sequence for the energy functional in V^nQ, 
which together with p.7[) and the fact that {$fc m .} converge to strongly in (L 2 (il)) N leads to 
(Aoo^oo) G 6oo, i.e., 

lim E($ k ) = £($00) = min E(V). 

Consequently, we obtain that each term of E($) converges and in particular 

lim ||V$fc m .|| ,n = HV^oollo.n- 
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Using p. 41) and the fact that Hq(SI) is a Hilbert space under norm ||V • ||o,n, we have 

Mm ||V($ fcm . -$00)110,0 = 0, 



which together with (|2.8[) . (|3.3p and (|3.5p implies (|3.6p . This completes the proof. □ 

To show that the limit in T^o n Q is indeed a ground state solution, we turn to the convergence 
of the a posteriori error estimators. We split the partition 71- into two sets T k + and Ti°, where 

T k + = {r G % : r G % V / > fc} and T k °=%\ T+. 

Actually, 71 + is the set of elements that are not refined any more, and consists of those elements 
that will eventually be refined. We denote by 

fi£ = U TeT +w fc (r) and fi£ = U TeT ouj k {T). 

Since the mesh size function h k = h k (x) associated with 7fe is monotonically decreasing and bounded 
from below by 0, we have that 

hoo{x) = lim h k (x) 

k — ^00 

is well-defined for almost all x G ft and hence defines a function in L°°(Q). Moreover, the con- 
vergence is uniform (see [13]), more precisely, if {/ifc}fc<=N is t ne sequence of mesh size functions 
generated by Algorithm 13. 11 then 

lim || h k - hoc || 0,00,0 = (3.8) 

k— > 00 

and 

lim ||/ifcXo?l|o,oo,o = 0, (3.9) 
where Xn° k i s the characteristic function of fl k . 

Lemma 3.2. Let (A^, G 0/j. If Assumption Al is satisfied, then there exists a constant C v > 
depending only on the mesh regularity, such that r]h(<&h,ft) < and 

Vh($h,T) < ||*fc||o,e,a; h (T) + ll^llwr) V T G 75,. 

Proof. Using (|2.15p , the inverse inequality, the Holder inequality, the trace inequality and Assump- 
tion Al, we have 

N N 

h T \\Tl T (^h)\\o,T = h T (22 II ~ ^ \j,h<l>j,h ~ 7,^4>i,h + V\oc<Pi,h + Vg\4>i,h 

N n 

< ^2 ^rdl^.hllo.r + ||A</>j,fc||o,T + || Vloc<^i,ft.||o,T + X! II II 0,r II <AU' 1 1 0,r 
i=l j=l 

+ l|eL(P*h)^»,/i||o,r + IK^" 1 * P<&J</>;,/Jo,t) 

< \\$h\\0,6, Uh (T) + \\$h\\l, Uh (T) 
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1/2 



and 

/ N \ 1/2 

fcV a ||J e (*/i)||o,e = ^ 2 (Ell2 V ^U^ + 2 V ^ h ^"^ ll ^J 

/ w 

< ^ /2 fE(H V ^l-llo,e + ll V ^UIIo,e) ] 

< ll*fc||l |Wh (T). 

Hence we obtain 

Vh($h,T) < ||$h||o,6,^(r) + ||*fc||l,o; fc (T) V T € Tft, 

which together with the Sobolev inequality implies i]h(&h,Q) < C r) , where the constant C v > 
depends only on the data and the mesh regularity. This completes the proof. □ 

Lemma 3.3. Let {$>k}keN be the sequence chosen by Alaorithm \3.1\ If Assumption Al is satisfied, 
then 

lim max ?yfe($fc,r) = 0. 

Proof. We see from Lemma \3. II that for any subsequence {3>fc m } of there exist a convergent 

subsequence {§k m .} and 4>oo € Ooo such that 

$fc m . -> $oo in W. (3.10) 

Now it is only necessary for us to prove that 

lim max rjk m . ($fc m ., t) = 0. 

i 

In order not to clutter the notations, we denote the subsequence {$fc m .}jeN by {$fe}fceN, and 
{7~k mj }jeti by {TfcjfcGN- We obtain from Lemma |3~21 that 

Vk(^k,T k ) < ||$fe|| ,6,^(r fc ) + ll*fc||x lW *(r fc ) 

< ||**-*oo||l i n + ||*oo||o,6 )Wfc (T 1k ) + ||^||l,ui»(7*), (3-11) 



where Tk € Mk be such that 



Vk($k,Tk) = max ?7fe($fc,T). 



Note that p. 101) implies that the first term on the right-hand side of (|3.11[) goes to zero. Since 
T k e Mk C T k °, we have from (j3"U)) that 

kfc(r fe )| < hl h < ||/ifeXngllo l0 o,n^° as fc ~> °°) 

which implies that the other two terms on the right-hand side of p. lip goes to zero, too. This 
completes the proof. □ 
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Define a global residual Rh(®h) € by 

AT W 

(Rh(3 fc ).r> = £ (ff» h ^ lh - ]T hj,h^,h,n) v r = (7,)^! e n. (3.12) 

We see that 

(R h (« A ).r)= E (fa-(*fc),r) T + E ( J «(*fc). r ) e ) vreH. (3.13) 

Thus 

|(R h ($ h ),r)| < £ 77h(*h,r)||r|| li(i;h(T) vren. (3.14) 

We shall show the weak convergence of H k (^ k ) in the following lemma. 

Lemma 3.4. Let {&k}keN be the sequence chosen in Algorithm ] 3. 11 If Assumption Al is satisfied, 
then 

lim (Rfc($ fe ),r) = VFeH. (3.15) 

fe— >-oo 

Proof. Similar to the proof of Lemma 13.31 it is only necessary for us to prove (j3.15|) for a extracted 
subsequence (still denoted by <£>&) that converge to <3>oo in H. 

Let T = (74)^! € (Hq(Q)) n and Ffc = (7i,fc)£Li € ^fc be such that 7^ is the Lagrange 
interpolation of 7,. We have from (|3.14l) that 

KR fc ($ fc ),r)| = KR fc (* fc ),r-r fc )| < £ %(**,T)||r-r fc || Wr) . 



For any e > 0, (|3.9I) implies that there exists n € N such that 

C n \\h n xnoJ\o,oo,n < £■ (3.16) 
Let k > n, we have C 7^, + C T k and hence 

|<Rfc($fc),r)| < Yl %($ fc ,T)||r-r fe || Wr) + £ »7fc(*k,r)||r-r fc || WT) 

rer+ rer k \T+ 

< r, k ($ k ,T+)\\r - r fc || lin + (T) + %($ fe ,r fc \ 7;+)||f - r fc || lin o (T) . 

Since Lemma yields 

%(^,r fe \r+) <%($ fe! r & ) < c„, 

we have from the standard interpolation estimate that 

KR fc ($ fc )»r)| <(vk($ k) T+) + cjh nX w n ||o,oo,n)||r|| 2 , n vre(ff 2 (fi)f, (3.17) 

Note that marking strategy (|2.18[) indicates 

rik($k,T+) < (#T+) 1/2 max % ($ fe ,r) < (#7;+) 1/2 max r), 
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which together with Lemma 13.31 implies that there exists N > n such that 

r) k (<Z> k ,T+) <e Vk>N. (3.18) 
By taking (13.161) . (|3.17[) . and (|3.18l) into accounts, we obtain 

lim (R fcm .($ fcm .),r) v r e (H*(n)) N , 

j— foo 3 3 

which together with the fact that H$(yi) is dense in H^fl) completes the proof of (|3.15p . □ 

Theorem 3.1. (convergence) Let {Qk\keti be the sequence generated by Algorithm \3.1\ If the 
initial mesh To is sufficiently fine and Assumption Al is satisfied, then 

lim E k = min£(*), (3.19) 
lim d n {Q k ,Q) = 0. (3.20) 

k— J-oo 

Proof. Let {(A^, *fc)}fcgN be the sequence chosen in Algorithm 13. II We know from Lemma I3TT1 that 
for any subsequence {(Afc m , *fc m )} mS N, there exists a convergent subsequence {(Ak m . ,&k m .)}j£N 
and (Aoo, *oo) <E Ooo such that 

*fc m . *oo in H, 
A km . ^Aoo in R NxN . 

3 

Consequently, it is only necessary for us to prove (A,*,,*,*,) £ 9, which implies (I3.19|) and (|3.20l) 
directly. For simplicity, we denote by {(A fe , *fc)}/c 6 N the convergent subsequence {{A km , , *fc OT . )}jeN, 
and by {Tk}keN the corresponding subsequence {Tk m 

We first show that the limiting eigenpair (A^, *oo) is also an eigenpair of (|2.6I) . We have from 
(|3~T2"j) that for any F £ H 

(H*«*oo-Aoo*oo,r) = (i?$ oo $ 00 -A 00 $ 00 ,r)-(R fc ($ fc ),r) + (R fc ($ fe ),F) 

= (i? $00 $oo - H^ k , F) - (A^^ - A fc $ fe , F) + (R fc (* fc ), T). (3.21) 

By a direct calculation using Assumption Al, we obtain 

(JT^Soo - < ll^ - *fc||x,n||r||i,n, 

which together with (|3 . 2 1[) leads to 

(fTft.^oo ~ Aoo^oo,^ < (11*00 - * fc ||i,n + lA^ - A fc |)||r||i in + (Rfe(* fc ),r). (3.22) 

Since Afe — > A^ and — > *oo in the first term on the right-hand side of (|3.22l) goes to zero 
when k goes to infinity. We obtain from Lemma 13.41 that the other term on the right-hand side of 
(I3.22|) goes to zero, and hence 

(ff* 00 * 00> r) = (A 00 # 00> r) wen. 

Then we shall show that for a sufficiently fine initial mesh, the limiting eigenpair (Aoo, *oo) is 
a ground state solution in 0. We set 

W = {(A,*) e M. NxN x H : (A,*) solves (l2~6l)l. 
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Note that 9 C W. Using the fact 



lim inf ||*-$|| in = VfeW, 



we can choose an initial mesh 7o such that 




Due to 7o C 71- , we have < Eq and hence (Aqq, $oo) € 8. This completes the proof. 



□ 



4 Quasi-optimality of adaptive finite element methods 

In this section we propose and analyze the following AFE algorithm using Dorfler's marking strategy. 
Algorithm 4.1. AFE algorithm with Dorfler Strategy 

1. Pick a given mesh To, and let k = 0. 

2. Solve (|2.16p on Tk to get discrete solutions ipi,k)(i = 1, • • • > AT). 

3. Compute local error indictors r/ki^k, t) for all r G Tk- 

4- Construct A4k C Tk by Dorfler Strategy and parameter 6. 

5. Refine Tk to get a new conforming mesh Tk+i ■ 

6. Let k = k + 1 and go to 2. 

We shall study the convergence rate and quasi-optimal complexity of Algorithm 14. 1[ for which 
we shall apply the perturbation arguments (c.f., e.g., [211 33 ) and certain relationships between 
nonlinear problem (|2.6I) and its associated linear boundary value problem (see (|A.ip ). 

To establish the relationship, we define 



a ($,r) = ^-(v0 4 ,v 7 «) v $ = r 



(70£i e n. 



One sees that there exists a constant c a > such that 



a(r,r)> Co ||r||? )n 



(4.1) 




a($,r) vre^, 



and K : W — > % be the inverse operator of £ such that 



o(A"#,r) = ($,r) vreE 

Note that (|4.1j) implies that K is well defined and there holds 

||ir*||i I n< ||$||-i,n Vfetf. 
Let Ph : % — > Vh be the -ff ^projection defined by 

a($-P h $,T) = V$6«,r€%. 



(4.2) 



(4.3) 



For any there hold 



||fl.*||i,n < ||*||i,n and lim ||$ - Ph^\\i,a = 0. 



(4.4) 
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4.1 Basic estimate 

First we recall an a priori error estimate, whose proof can be referred to |13) . Define 

x^ h = s NxN x (v h n (5» © 7i)). 

Theorem 4.1. Lei (A, $) 6e a solution of (|2.6I) . If Assumptions A2 and A3 are satisfied, then 
there exists 6 > smc/j that for sufficiently small h, (|2 . 14[) /ias a unique local solution (A/,,$h) G 
AT$./i H -B,5((A, $)). Moreover, we have 

||4-fcfc||i,n< mf ||$-*||i,n, (4.5) 

|A/, - A| < p fc - $||?, n + - *|| ,n, (4.6) 

||*-*h||o,n <r(/i)||*-* h ||i, n (4.7) 

luii/i r(/i) — > as /i — > 0. 

Using Theorem 14. 11 we can denote afterwards by (A^, $^) € X^/i n Bj((A, <&)) the unique local 

discrete approximation of (A, $) e 9. We say an equivalence class \$>h] approximate the equivalence 

class [$] if there exists an orthogonal matrix Uh such that $^ = ^>hUh- For simplicity of notation, 

t ( \ 

we denote by V = V\ oc + V n \ and N(p<s>) — / i ;dy + e ' xc {p^)- 

Jn\'-y\ 

Lemma 4.1. Let (A, $) be a solution of (|2.6p and ho £ (0, 1) be the mesh size of the initial mesh 
Tq. If Assumptions A2 and A3 are satisfied, then there exists k(h) such that k(h) — > as h — > 
and 

\\V($ h - *)||_i, n + ||A/"(p*J$h - A/"(p*)$||_i,n < - $ h ||i,n. (4.8) 

Proof. For any 'J £E "H, by using the Holder inequality and the Young inequality, we have that for 
any e > 0, there holds 

||*||o,3,o < M^llo^Il^llS^^ = C^- 2/3 M*llo^C^ 2/3 M^llg^«) 

< -3-11*110,0 + ^11*111,0, 

which together with (|4.7p implies that there exists a positive constant C independent of h and e 
such that 

||*-*A||o I 3,n< (Ce- 2 r{h)+e) ||$-$ h ||i,n Vft6(0,M- 
Taking e = r(Zi) 1 / 3 and setting = r(Zi) 1 / 3 , we have that k(h) — > as h — > and 

||$ - $fc||o,3,o <k(/i)||*- $/,||i,n, 
which together with the Holder inequality leads to 

||Moc($ - *fc)||-i,n = sup &^LZ^l£) < ||Moc||o,o||$ - *fc||o,3,n < fc(>OII* - (4.9) 
rew II 1 ||i,o 

For the nonlocal pseudopotential operator, we derive 

(Ki($/, - r) < ||$ - $h||o,n||r|| ,n vrew 
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from the fact that 



(X](CiAh < ||&,fc-&||o,nNlo,n V^e^fl), i = l,--- ,N. 

3=1 

Therefore, we have 

\\V nl ($ h - $)\\-i,n = sup (M^Jltlll) < ||$ _ $ fe || o n < r ( fc )||$ _ (4.10) 
rew II 1 II i,o 

For the exchange-correlation part, we have that there exists £ = - ■ ■ ,£jv) G "H with £j = 
$i4>i,h + (1 — ^i)^ an d <5< € [0, 1] (i = 1, • • • , iV), such that 

N 

«c(p*J*fc ~ e^ c (/9*)$,r) = ^ / (4c(P?) + 2^ 2 exc(P?))(^,ft " &)7v 
This together with Assumption A2 leads to 

JV 

(e' xc (pz h )$ h -e' xc (pz)$,r) < J2 / {Pt+Pf)\<l>i,h-<t>i\-\li\ 

i=i Jf2 

AT 

~ X! (ll/ , f t llo,3/a,nll<Ai,h - ^i||o,n||7»llo,6/(3-2o),n + ||Pillo,3,n||0i,h - ^i||o,n||7i||o,6,n) 

i=l 

< ||*fc-*||o,n||r||i, n VTeW, (4.11) 
where the Holder inequality and the fact 

l|P€llo,3,n < ||Cl|g,e,n < ll«ll§,e,« -h ll*fcllg,«,n < C- 

are used. For the Coulomb potential, we obtain from the Young's inequality and the Uncertainty 
Principle 05] that 

JV 



\\r 1 * 0$ - p*J||o l0 o,n < XI H v (^ + &Vi)llo,£i||& - <M|o,n < II* - *h||o,n- 

i=l 

Therefore, we have that for any v € -£?o(^) and 1 < i < N, there holds 



<2 



< ||r _1 * p$J|o,«>,n|l0i,/> - 0i||o,n|M|o,n + Ik -1 * (p*,, - />*)||o,oo,n||0i||o,n|Mlo,n 

< ll^i - ^i,/j||o,n||f ||o,n + II* - *h||o,n|Mlo,fi) 
which implies 

((r-^PttJiA-fr- 1 */»»)$, T)< ||* -$ h || ,n||r||o,n V r e W. (4.12) 
Consequently, we obtain from (|4.1ip , (|4.12[) and the definition of N that 

Wip*>)*H-Mif)*U* = *9 (^)^~^) $ ' r ) < ||* _ * fc | knj (4.13) 

reH II 1 ||i,n 

Taking into accounts (14771) . (|49]) . (|4~T0l) and (l4~T3|) . we complete the proof of 1(378]) . □ 
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We now exploit the relationship between the nonlinear eigenvalue problem and its associated 
linear boundary value problem, which will be employed in our analysis. We rewrite (|2.6[) and (|2.14l) 

as 

$ = K ($A - V$ - 7V(p<t,)$), 

$ h = P h K($ h A h - V$ h - N( P ^ h (4.14) 
respectively. Set W h = K($ h A h - V$ h - Af(p<s> h )$ h ), we have $ h = P h W h . 

Theorem 4.2. Let (A, $) oe a solution of (|2.6p . If Assumptions A2 and A3 are satisfied, then 
there exists n(h) € (0, 1) suc/i i/iai — > as /i — > and 

||*-*fc||i,n = \\W h -P h W% <n + 0{K{h))\\$-$ h \\ 1 ,a. (4.15) 

Proof. By the definition of W h , we have 

$ - W h = X($A - $ h Afc) + KV($ h - 4) + K(AA(p*J$ h - AA(p$)$). (4.16) 

For the first term on the right-hand side of (|4. 1 6[) . we obtain from (|4.2I) and (|4.7I) that 

\\K($A-$ h A h )\\ lin < ||*A - $fcA h ||o,n < ||(*- */.)A|| ,n + ||*h(A - A fc )||o,n 

< ||*-#h||o,a|A| + |A-A h |<r(fc)||*-$h||i,n. (4.17) 

Using Lemma [4~T1 we can estimate the second term on the right-hand side of (|4.16[) as follows 

||ifV($-4 A )||i jt i < ||V(* — * A )||_i, n < /cWH*-*fc||i,n. (4.18) 
Using (|4~2]) . ([4~7l> and (j478|) . we obtain for the last term of (|4~T6|) that 

||if(A/'(p*J*fc-A/'(p«)$)||i 1 n < ||A/'(p«J$fc-A/'(p«)$||_i 1 n<r(A)||4-*fc||i,n. 
Set K(h) = r(/i) + k(/i), we derive from (l4~TtJj) . (|4TT7j) . (|4TT51) and (|4TT5j) that 

||$ - ^Hx.n < - $ h ||i, n , (4.19) 

with C* being some constant. Note that (14. 14[) implies 

$ - $ h = W h - P h W h + $ - W h , 
which together with (|4.19p leads to (|4. 1 5|) . This completes the proof. □ 

4.2 A posteriori error estimates 

Define 

k(h ) = sup n(h) (4.20) 

he(Q,h ] 

and note that h(ho) <C 1 if ho <C 1. Based on the relevant results for linear boundary value problems 
(see Appendix), we have the following a posteriori error estimates for AFE approximations. 
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Theorem 4.3. (a posteriori error estimate) Let (A, <£>) be a solution of (|2.6[) , ho <C 1 and ft, 6 
(0, ho]. If Assumptions A2 and A3 are satisfied, then there exist positive constants C\,Ci and C 3 
depending on the coercivity constant c a ( in (|4.1|) ) and the shape regularity constant 7* ( in (|2.12j) ), 
such that 

II* -*h||? j0 < CmU^n), (4.21) 

C 2 ^(* h ,n) < ||$-$ ft ||? i0 + C 3 os4($ h ,fi). (4.22) 
Proof. Due to £W / ' 1 = $ h A h - - J\[{p<s, h )<$> h , we obtain from (|A.7|) and (|A.8j) that 

||W* - P h W h \\l n < C 1 f h (P h W h ,il), (4.23) 

^(P/J^O) < \\W h - JW h ||?,n + C 3 5&l(P h W h ,n), (4.24) 

where the constants C\, C2 and C3 are given in Theorem lA.il fj^{PhW h , fi) and 6sc^(P/ l W / ' 1 , f2) are 
defined by (|A.5[) and (IA.6I) with T being replaced by Py t W h . It is easy to see that rjh(PhW h , fl) = 
r/ h ((S> h ,n) and 6sCh(PhW h , 17) = osc/^/j, Q) from their definitions and the fact that & h — PhW h . 
We have from (j4~l%)) and (|4T20|) that 

11$ - *fc||i,n < (1 + Ck(h Q ))\\W h - fWlkn, 
which together with (I4.23[) leads to (|4.21l) by taking the constant 

d =C 1 (1 + Ck(h )) 2 . (4.25) 
Similarly, we get (|4.22l) from (|4.14l) . (|4.15l) and (|4.24|l . In particular, we may choose C 2 and C 3 

by 

C 2 =C 2 (l-Ck(h )) 2 , C 3 =C 3 (l-Ck(h )) 2 . (4.26) 
This completes the proof. □ 

Our analysis is based on the following crucial technical result. 
Lemma 4.2. For given constant 8 £ (0, 1), if 

E vt(*k,T)>6 V t(* k ,n), (4.27) 

T£M k 

then for any $>k = ^kUk with Uk being some orthogonal matrix, there exists a constant 9' £ (0, 1), 
such that 

E 4{*k,T)>0'rrl(* k ,n). (4.28) 



T£M k 



Proof. We write Uk = ( a ij)- Since [4 is orthogonal, we have J2f=i { a i,jY = ^ ^ or any * = 1> ' " ' > ^> 

and Eili ( a ij) 2 = 1 for any j = 1, • • • , iV. 

On the one hand, we obtain from = ^kUk that 
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Consequently, we have 

N N N N 



4^ k , n) = J2 ^(E < E ( E a "M^,k, (i)] 

i—l j — 1 2—1 i— 1 

N N N N N N 

^ E(EK*) 2 E^^)) =EEfc' fi )^Efc> fi )> 



where the fact X^=i( a £i) 2 = 1 is used. Namely, 

r,l($ k ,Sl) < Nr,l(* k ,Sl). (4.29) 
On the other hand, $fe = ^kUk implies ty k = ^kU^ . Hence, 

JV 

ipi,k = E a u^'" i = l,---.^V- 

By the similar process and the fact Yl,j=i{ a i j) 2 = 1 we obtain that 

E vlm,r)<N Vt(^k,r). (4.30) 

r£M k reM k 

Combining (|4~27|) . (|4~29| and (f4~30)) . we arrive at 

4{9 k ,n)<N4(^ k ,n)<N^ ^ 4(* k ,r) < n 2 ± J2 ^< r )< 

which is nothing but (|4~28| with 9' = □ 
4.3 Convergence rate 

We turn to analyze Algorithm 14.11 We shall first establish some relationships between two level 
finite element approximations. We use Th to denote a coarse mesh and 7~h to denote a refined mesh 
of T H - 

Lemma 4.3. Let h,H E (0, ho] and (A, $) be a solution of (|2.6|) . // Assumptions A2 and A3 are 
satisfied, then 

\\$-$h\\i,n = \\W H - P h W H \\ 1 ,n + O(k(h ))(\\<P-<i>H\\ h n + ||$ - *h||i,n) , (4.31) 
osc h ($ h ,n) = 6Sc /l (P /l M/ ff ,fi)+0(K(^o))(||$-^||i,o + ||$-$h||i,n), (4.32) 

and 

= %(P/ l V^ ff ,r>) + O(K(^ ))(||$-$H|| 1 ,o + ||$-$h||i,n). (4.33) 
Proo/. First, we obtain f)4.31[) from (|4~4|) . (|4.19|) and the identity 

$ - $ h = VF H - PftW^ + Pfc(W H - W h ) + $ - VF ff . 
For the estimate of ([432")) . we get from = P ;i VF ff + P h {W h - VF ff ) that 

oTc^P^, SI) < 5Sb h (P h W H , SI) + 5Sc h {P h {W h - W H ), SI), (4.34) 
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where osc is given in Appendix. Using (I4.14[) and the fact osc? l ($/ t , 0) = osc/,. ($?,,, ^), we know that 
it is only necessary to estimate 6sCh{Ph(W h — W ), O). 

Since ^CPU' 1 = $ fe A ft - V$ h -Af(p<s> h )$ h and £VK H = $ H A H - - Af(p$ H )§h, we obtain 

- W H ) = <Z> h A h - <S> H A H + V(<S> H - $h) + Ar(p* H )$H -^(p»J$fc. 

Let G = P^fW 1 - JU H ) and ft T (G) be defined by (TOt with T being replaced by G. We have 

£r(G) = $ h A h - <S> H A H + V(<S> H - +Ar(p* H )*H-tf(j>* h )*h - £G 

and 



5Sc 2 h (P h (W h -W H ),il) = ^(G,r) = ]T ^^(G) -^r(G)||g, 



< J2 K\\Kt(G) + CG- (K T (G) + CG)\\l T + J2 %\\£G - CGf 0>T . (4.35) 

Using the inverse inequality, and the fact that $hAh and §hAh are piecewise polynomials 
vectors over Th and Th respectively, (|4.7[) . and f|4.8[) . we have 



.1/2 



( J] tf||£r(G) + £G - (£ T (G) + £G)||2 T ) 

< hT(\\V($ H -$ h )\\ ,r + \Mp* H )$H-M(p* h )$h\\0,T) 

< ft(/»o)(||*-*h||i,n + ||*-$ff||i,n). (4.36) 
Combining the inverse inequality, (|4.4p and (|4.19[) . we arrive at 

( Y, K\\CG-CG\\l T ) 1/2 < ( 2 ^II^IU V2 S »,n 

= \\Ph(W h -W ff )||i,o<K(/io)(||$-^||i,n + ||$-$H||i,n). (4.37) 

Taking ((4735]) . (|4736| and (j4~37| into account, we obtain 

ffic^P,^ - W H ),Q) < k(Ho) (||* - & h \\i,n + II* - *ff|M , (4.38) 

which together with P~3"Ij) leads to (|4~52l . 

Finally, we shall prove (|4~3"3"|) . We obtain from (|7Q|) . (gUg]) and (|273"gj) that 

Vh{Ph(W h - W H ),n) < \\(W h -W H )- P h {W h -W H )\\ hn +5§c h (P h {W h ~w H ),n) 
< k(h ) (||$ - #fc||i,n + II* - *ff||i,n) • 

This together with the fact 

^(iw' 1 , n) = fi h {p h w H + p i (VK /l - vu ff ), n) 

leads to 

%(PhW h , fl) = %(PfcW ff , fi) + 0(£(M) (II* - *ft||i,n + ||* - *H||i,n) , 
which is nothing but (|4.33p . This completes the proof. □ 
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Applying Theorem IA.2I and Lemma T4. 31 we get the following theorem. 

Theorem 4.4. (error reduction) Let 9 6 (0,1) and ho <C 1. Let {^k}kefi<, be a sequence of finite 
element solutions corresponding to a sequence of nested finite element spaces {Vk}k<£N produced by 
Algorithm ^. 1\ Assume [$fcj is an approximation of some [$] with $ being one solution of i2.6\) , 
denote ki + \{> ki) the minimal index among all indexes fc(> ki) which satisfy that \^k\ approximates 
[<&]. If Assumption A2 is true and (A, $) satisfies Assumption A3, then 

II* - ®k i+1 \\U + 7Vk i+1 (®k i+1 ,T ki+1 ) < e (||$ - **J?,n + 7»Jk 4 (**,>7i t )) (4.39) 

with <&ki € [*fcj H A$ : fc i; 7 > and £ € (0,1) some constants depending only on the coercivity 
constant c a , the shape regularity constant 7* and the marking parameter 9. 

Proof. Recall that $fe i+1 G [^a^J H A$ : fc i+1 and e [^fcj n -X*,^ satisfy the a priori error 
estimates (|4.5p and (|4.7p when /i is replaced by and fej. For convenience, we use to 
denote $fc i+1 and $fci, respectively. Then it is sufficient to prove that for <&h and <f>#, there holds, 

II* - *&lli,n + 7^(*fc,n) < £ 2 (ll* - *^lli,n + 7^(*H,n))- 

We see from Lemma I4T21 that Dorfler Marking strategy in Algorithm 14.11 implies that there exists 
a constant 9' — € (0, 1), such that 



2^ ^($ ff ,r)> 0^(^,0). 



Therefore, we obtain from W H = K($ H A H - - N(p9n )*ff) and $ H = P H W H that Dorfler 
strategy is satisfied for W H with = -Jy . So we conclude from Theorem IA.2I that there exist 
constants 7 > and £ € (0, 1) satisfying 

\\W H - P h W H \\l a + jf h (P h W H , fi) < f 2 (||T4^ - $ H ||; it , + 77&($* , ( 4 - 4 °) 

where the fact f)n{PnW H , fi) = r] H {^H, SI) is used. 

From (|4.19|) . we get that there exists constant Ci > such that 

'l + CiR(/»o)) ||$ - *H||;, n + 7»&(*tf> O) > ll^ H - PffW ff ||?, n + 7^1(^,0). (4.41) 

By Lemma 14.31 and the Young inequality, we have that there exist constant C2 > such that 



II* - *fc||i,n + lvU$h,Sl) < (1 + h)\\W H - PhW H \\l n + (1 + 5 1 )jf h (P h W B , SI) 

+C 2 (1 + S^)K 2 (h )(\\^ - <S> h \\l n + ||$ - * ff ||?,a)- (4-42) 

where o\ € (0, 1) satisfies (1 + <5i)£ < 1. 

Combining ([4~4TJ]l . (l4T4Tp with g^U, we get that 

(l - tf 3 (l + JrVCk))) II* - */i||?, n + T7 2 (*fc, O) 

< ((i+fc)f a +(i+*i)^6i/KM+&(i+^ 



Since /i -C 1 implies k(ho) -C 1, there holds 

J7 

l-Cs^^^o) 



* - *hH 2 o + ; a r -i- 2 „ : ^(^» fi ) 



1 - Ca^-^ho) V (l + ^ + C^M 
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with C3 some constant depending on C\ and C%- Note that ho -C 1 implies k(ho) < 1, we see that 
the constant £ defined by 



(1 + s^e + c 3 k(h Q ) 



1/2 



satisfies £ € (0, 1) when /i « 1. 

Finally, we arrive at (I4.39|) by using the fact that 



<7 with 7 = » 2i. ■ ( 4 - 43 ) 



(1 + 5i)| 2 + CaftCM l-C^f^M' 
This completes the proof. □ 

We have from Theorem 13.11 that if {^k} is obtained by Algorithm 14. 1[ then there exists a 
subsequence {[^fcj} that converge to some equivalent class [$], where $ is a solution of (|2.6[) . 
Here, a sequence {[^fcj} converges to a equivalent class [$] means that there exist orthogonal 
matrices J7 fci £<D NxN , such that 

lim V ki U ki = 

z— f 00 

Therefore, combining Theorem I4.4| we can easily have the following theorem. 

Theorem 4.5. (convergence rate) Let 8 £ (0, 1) and ho *C 1. £e£ {^fc}fc e N ^ e a sequence of finite 
element approximations obtained by Algorithm [^77] and {[^fej} 6e </ie subsequence that converges 
to some [$], where $ is a solution of K2.0i) . If Assumptions A2 and A3 are satisfied, then there 
exists e [^fej] H ^"$.fei /or aH i € N, suc/i £/ia£ 

II* - ^ i+1 ||i,n +-yvl +1 ($k i+1 ,T ki+1 ) < e (||0> - $ fci ||ta + 7<(^ iJ Tfe < )) , (4-44) 

where 7 > and £ G (0, 1) are constants depending only on the coercivity constant c a , the shape 
regularity constant 7* and the marking parameter 8. 

4.4 Complexity 

Finally, we study the complexity of Algorithm 4.1 in a class of functions. Define 

^={feH: |*| s , 7 < 00}, 
where 7 > is some constant and 

|*k 7 = su P£ mf (#r-#r ) s 

£>0 {TCTo: int4, re v T (||*-*rll?, n + (7+l)os4(* r: r)) 1 /2< e } 

and T C To means T is a refinement of To- It is seen from the definition that, for all 7 > 0, 
= A{. For simplicity, here and hereafter, we use A s to stand for A\, and use |*| s to denote 
|^| a,7' So A s is the class of functions that can be approximated within a given tolerance e by 
continuous piecewise polynomial functions over a partition Tk with number of degrees of freedom 
satisfying #T fe - #T < e-V'\*$ la '. 

Lemma 4.4. Suppose 8 £ (0, 1) and ho <C 1. £e£ and $ft, &e i/ie solutions of ()2.16j) over a 
conforming mesh Th and its refinement Th, and h] an d [&h] approximate the same solution class 
[<f>], where $ is a solution of H2.6]) . Suppose that for some <j> € [<&] satisfying (|2.11|) , we Ziaue 

II* ~ <Ma,n + 7*os C 2($,, n) < #(||$ _ + 7 ,osc|,(<i>H, fi)) (4.45) 
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with $ h € n X$ ih and $ H € n X$, H , 7* > and e (0, TTien, i/ie set 

7?. = IZth^Th satisfies the following inequality 

ren t£T h 
here 8 = A fA^)!"^'!; s- N j wii/i Cq,0*, and 7* &emn constants defined in the proof. 

t'o(t'i + (l+' i W<-'i)7*) 



Proof. For W H = K[$ H A H - V$h - N(p® H )®h) , we observe from Lemma H31 that 

11$ - * h ||i,n = \\W H - P h W H \\ hQ + O(k(h )) (\\W H - P„W H \\ hn + \\W H - P h W H \\ hn ) , 

osc h ($fc, f>) = 6sc h (i\W H ,fi) + 0(k(/»o)) - PffW^llLn + \\W H - P H W H \\i, a ) , 

Proceeding the similar procedure as in the proof of Theorem 14. 4[ we have 
\\W H - P h W H \\l n +%6§-c 2 h (P h W H ,n) < fc(\\W H - P H W H f a< n + %^ 2 H (P H W H ,n)) (4.46) 
with 

\ 1 - c 4 ar K 2 (h ) ) 1 - c^r K 2 (h y 

where C4 is some positive constant and 6\ G (0, 1) is some constant as shown in the proof of 
Theorem 14.41 

Set Co = max{l, Sa}, we get from (|A.8|) that 

(1 - 20, 2 )C 2 f 1 2 H (P H W H , fi) < (1 - 2/3?) - PirW^ l^n + C 3 6Tc^(Ph^, O)) 

< C (l - 2/3?) - iW*||2,n + %6s&(P*W* fi)) , 

which together with (I4.46[) produces 

^(1-2/3?) £ i&CiWV) < (\\W H -P H W H \\l n +%5^ H (P H W H ,n) 

- \\W H -P h W H f a>n - 2%oscl(P h W H ,n)). (4.48) 

Thus using equality 

\\W H - P H W H \\l n -\\W H - P h W H \\l n = \\P H W H - P h W H \\ln 



and Theorem IA.3[ we obtain that 

\\W H - P H W H \\l n -\\W H - P h W H \\l n < Cx £ Vh(PhW h 7 t). (4.49) 

By the triangle inequality, the inverse inequality, and the Young inequality, we can easily see that 

]T 5Sc 2 H (P H W H ,r)<2 Yl 5Sc 2 h (P h W H ,T) + 2C 2 jP„W H -P h W H \\l n 
T£T H nT h T£T H nT h 
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and 

5Sc 2 h (P h W h ,t) <fj 2 H (P H W H ,T) VreTff, 

where C* is a positive constant depending on the shape regularity constant 7*. 
Hence, we may estimate as follows 

6s<? H {P H W H , Q) - 2osc h (P h W H , Q) 

Ten T£T H nT h T£T H nT h 

< E V 2 h(PhW h ,t) + 2C 2 \\P H W H - P h W H \\l n 

< (1 + 2C 2 C 1 )J2vh(PhW h ,t). (4.50) 

Ten 

Combining (|4.48[) . (|4.49|) and (|4.50[) . we then arrive at 
namely, 

E^(^> t )>^ E v%($h,t) 

Ten TeT H 

with 

* C 2 (l - 2/3?) 



C (Ci + (l + 2C?Ci)7,) 

This completes the proof. □ 

We mention that Dorfier Marking Strategy selects the marked set Aik with minimal cardinality. 

Lemma 4.5. Let 8 G (0, 1) and ho <C 1, {^fcjfceNo be a sequence of finite element solutions 
corresponding to a sequence of nested finite element spaces {Vk}keN produced by Algorithm \4-l\ 
and 6 G (0, c3{C\+(\+2C' i Ci)i) )■ If [^fc] approximates the solution class [$], where $ is a solution 
of H2.6\) , then for any $ G [$] H -4 s satisfying (|2.11[) . we /iaue 

< (||$ - n + 7o^($ fe , fi))~ 1/2s |$|y s , (4.51) 

where $fc G [\I/fc] H-X^fc, and £/ie hidden constant depends on the discrepancy between the marking 

parameter C3(Cl +a+2C^c 1 ) 7 ) and 61 ■ 

Proof. Let a, ai G (0, 1) satisfy ai G (0,a) and 

< C 3 (C 1 + (1 + 2C2C 1 )7) 1 " J ' 
Choose <5i G (0, 1) to satisfy (1 + <5i)£ 2 < 1 and 

(1 + 5ifa\ < a 2 , (4.52) 
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which implies 

(l + <5i)af<l. (4.53) 

Set 

e = - *k\\l,n + Josct(<P k ,n)) 1/2 

and let T E be a refinement of 7o with minimal degrees of freedom satisfying 

ll*-*«ll2,n + (7+l)°"£(* e ,n) <e 2 . (4.54) 

We get from 3> G *4 S that 

#Te-#T < e - 1/s |<&|i /s . 

Let 7Z be the smallest common refinement of Tk and 7^. Note that W £ = K(§ £ A £ — V$ £ — 
N{p<s> € )$ e ), we obtain from the triangle inequality, the inverse inequality, and the Young inequality 
that 

oa£{P*W e ,tl) < 25S-cl(P £ W e ,n) + 2C 2 jP £ W s - P*W E \\l n , 
where P £ and P* are Galerkin projections on T e and 7* defined by (|4.3[) . Due to identity 
ll^ £ - P^lllo = ||W e - P e W% n - \\P*W S - P e W%n, 

we have 

W e -P*W e \\l u + ^™cl{P*W e ,n) < \\W s -P s W% Q + -^osc 2 (P e W £ ,Q). 
Since ()A.9|) implies 7 < -^jz-, we obtain that 

||W £ -PJ^ ;n +7^fc*(P*W £ ,fi) < \\W s -P e W s \\l u + -^osc 2 s (P s W e ^) 

< \\W £ - P £ W s \\l n + (7 + a)o SC 2(P £ ^ e , n), 

where <J — — 7 € (0,1). Applying the similar argument as that in the proof of Theorem 14.41 
when (14.33[) is replaced by (|4.32[) . we arrive at 

||#-**|£n + 7°«2(4.,fi) < a 2 (\\$-<S> £ \\l n + (j + <j)o S c 2 £ (P £ W^n)) 

< al(\\$-$ £ \\l n + ( 7 +l)osc 2 (P £ W e ,n)), (4.55) 

where 

a2 = (1 + Sj) + C 3 k(h a ) 
a ° l-C 3 S^k 2 (h ) 

and C3 is the constant appearing in the proof of Theorem 14.41 It then follows from (|4.54p and 
(|4~55|) that 

II* - **\\lfl+1os<2($ m ,%) < a 2 {\\$ - $ k \\la +loscl{® k ,%)) 
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with a = A=a ai. Using (14.531) . we obtain a 2 € (0, i) when h -C 1. Set 9 = - ,/^ 2 ^ , ^° 2 a NM , 

7 = , \ ' = max(l, %), and d 2 = (1 1 + % ) f-t^» ( ^" ) ■ We have from Lemma EI that % 

satisfies 

where 7Z = 7Zj- k ^j- r is the refined elements from 7fc to T*. 

It follows from the definition of 7 (see (|4.43|1 ) and 7 (see (|A.9[1 ) that 7 < 1 and hence Co = 
Since h < 1, we get that 7 > 7 and a £ (0, ^=a) from (|3~52"]) , We see from (|4~2"5j) . P~2f))) and 
7 > 7 that 

^(i^*) > , . ^ _(W) 



f (^ + (1 + 2(72(71)7) (73(^ + 1 + 2(72(7!) 

c 2 



(1-C£(M) 2 , _ 2) 



(l-CB(ft )) 2V i((l+Cs(ft )) 2 ) * (1 + Ck(h )) 2 ' 

- C7 3 (f + (1 + 2C2d)) (1 ~ = ^(d + (1 + 2C?C 1 ) 7 ) (1 " ^ > 



when /iq -C 1. Thus we arrive at 



#M k < m < - #r fe < #7; - #r 



< (^«i)^ 1/s (II* - <Ma,^ + 7 OS c 2 ($ fe ,r fe ))^ 1/2s i$iy s , 



which is the desired estimate (14.51)) with an explicit dependence on the discrepancy between 9 and 
C27 

C* 3 (Ci + (l+2CJCi) 7 ) 



via ax- This completes the proof. □ 



Theorem 4.6. (optimal complexity) Let g (0, 1) and ho <C 1. Assume that Assumption A2 is 
satisfied and not accounting the invariance of unitary transform (|2.6I) /ias to solutions, which are 
denoted as [<f>"'](l = 1, • • • , to) where to can oe chosen to be 00. Let {^fcj-fceNo ^ e a sequence of finite 
element solutions corresponding to a sequence of nested finite element spaces {Vk}kefS produced by 
Algorithm ^. 1\ Then the following quasi-optimal bound is valid 

m — x/2s 

#r„-#To<^(||cl> i -$^|| 2 !0 +70sc 2 iii ($^,r2)) , (4.56) 
1=1 

where € [$^]n.4 s satisfying (|2.11j) . <& l k e iH-Xgi^ an d the hidden constant depends on 
the exact solution $' and the discrepancy between 9 and c 3 (Ci+{i+2C 2 Ci)f) ' Here, n.; and fc ni are 
</ie toiaZ number and the maximal index of iteration which approximate [$^](Z = 1, • • • ,to) among 
the n iteration, respectively. 

Proof. Assume that among the iterate solution spaces {^i}" = i, there are n\ approximations for 
[$(')](/ = 1, • • • , to), which are denoted by ^P/^, « = 1, • • • ,n\. Here, X)I=i n i ~ n > ano - n i can be 0- 
Recall that (see Theorem 6.1 in [54 ) 

m rti 

#T n -#T <J2I2# Mk <> 

1=1 1=1 
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we obtain from (|4.51p that 

rn m 

#r„ - #r < EE - <M,n + tosc^^,^))- 1725 (\&\i /s )- 

1=1 i=l 

Note that P~2"2"]l implies 

II*' - < II l,n + Til « ,n)<C(\\&- $i Win + 70sc 2 fei (^ , fi)) , 
where C = max(l + ^-). We then arrive at 

m m 

#T n -Wo < EE(di* , -*i.ii?.n + 7»jg l (*i„n)))" 1/a '(i*'i; / '). 

1 = 1 i=l 



Due to (|4.44|) , we have that 

II*' - <, ||?,n + 7< «, . < e 2( " ! - 4) (11$' - $L II i,n + 7^ , «)) 
Consequently, we obtain 



< E(i $ ^ /s (ii^- $ Lji^+7<(C,o))- 1/2s Er 



< 



E (is'i^Oi* 1 - *L, n 2 ,n + 7< «, , n)) 



-1/2; 



where the fact £ < 1 is used 

?i,n)<Tfr(4 fc 



Since osc k ($ l kl Q) < 7? fc (<I>j., D), we get 



m — l/2s 

#7;, - #T < E (V - *L, H?,n + 70"!, , n)) 
z=l 

This completes the proof. □ 

5 Numerical examples 

In this section, we shall present some numerical simulations for three typical molecular systems: 
CgHgOi (Aspirin), Cs-ffgC^A^a amino acid), and 6*60 (fullerene) , which support our theory. Due to 
the length limitation for the paper, we only show the results for pseudopotential approximations 
for illustration. 

Our numerical experiments are carried out on LSSC-III in the State Key Laboratory of Scientific 
and Engineering Computing, Chinese Academy of Sciences, and our package RealSPACES (Real 
Space Parallel Adaptive Calculation of Electronic Structure) that are based on the toolbox PHG 
[66] of the State Key Laboratory of Scientific and Engineering Computing, Chinese Academy of 
Sciences. 

In our computations, we use the norm-conserving pseudopotential obtained by fhi98PP soft- 
ware and the LDA exchange-correlation potential. We use Algorithm 4.1 and apply the standard 
quadratic finite element discretizations. Since the analytic solutions are not known even for the 
simplest systems, we only show the convergence curve of the a posteriori error estimator ^(^fc, fi) 
in our figures. The mesh and density illustrations are drawn using ParaView. 
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Example 1: Aspirin CqH^O^. 

The ground state energy obtained by SIESTA is —119.621 a.u.. In our computations, we choose 
the computational domain to be ft = [—20.0, 20.0] 3 . 

The atomic configuration, the calculated ground state charge density and the associated com- 
putational mesh are shown in Figure 15.11 First, comparing the configuration figure (the left one 
of Figure ISTTl) and the charge density figure (the middle one of Figure IBTTj) . we can see qualita- 
tively that our calculations are correct, the carbon-hydrogen bonds, carbon-oxygen bonds, and the 
oxygen-hydrogen bonds are preserved very well in our calculations. If we take a detailed look at the 
charge density figure, we can further see that the charge is more concentrative around the oxygen 
than around the carbon. We also see from the mesh figure (the right one of Figure 15. 1|) and the 
charge density figure that our error estimator can catch the oscillations of the charge density very 
well, which qualitatively implies that our error estimator is efficient. 





Figure 5.1: CgHgO^: configuration, charge density and mesh on plane z = 0. 



We now turn to analyze some quantitative behavior of our calculations. The convergence curve 
of the ground state energy is shown in the left of Figure 15.21 We observe that the ground state 
energy approximations converge to —119.918 a.u., which is very close to the value given by SIESTA. 
This result validates our calculations quantitatively. We see from the right of Figure 15.21 that the 
convergence curve of the a posteriori error estimator is parallel to the line with slope — | , which 
means that it reaches the optimal convergence rate. From the analysis result for the a posteriori error 
estimator (Theorem 4.3) the optimal convergence of the a posteriori error estimator also indicates 
that the approximation of the eigenfunction space have reached the optimal convergence rate, which 
coincides with our theory in Section 4. 



10000 100000 le+00 

degrees ol freedom 



10000 100000 le+00 
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(a) ground state energy (b) r/h^t^Q) 

Figure 5.2: The convergence curves of the ground state energy and rj^i^h, 
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Example 2: a amino acid Cs-ffgC^iV '. 

The ground state energy obtained by SIESTA is —75.494 a.u.. In our computations, we choose 
the computational domain to be fl — [—10.0, 10. 0] 3 . 

The atomic configuration, the calculated ground state charge density and the associated compu- 
tational mesh are shown in Figure [5731 We have to point out that for Cs-ffgC^A , not more than 2 
atoms stay in the same plane. Therefore, it is very difficult to find a plane where the configuration 
and the charge density coincide very well with each other as Example 1. Similar to Example 1, 
we also choose the plane z = as our viewpoint. Anyway, we can see from the figure for charge 
density and the figure for the adaptive mesh that our error indicator is very efficient. These results 
can validate our computations. 




Figure 5.3: C5HQO2N: configuration, charge density and mesh on plane z = 0. 



The convergence curves of the ground state energy and the a posteriori error estimator rjk^ki ^0 
obtained by the quadratic finite elements are shown in Figure 15.41 from which we observe that the 
ground state energy approximations converge to —75.494 a.u., and the a posteriori error estimator 
decays with a rate — |. This implies the similar conclusions as those for Example 1. 



-80 I ' ' ^ ' ' -J ' ' ^ ' ' ^ ' ' ^ 0.1 I ' ' — ' ' — ' ' — ' ' — 

100 1000 10000 100000 1e+06 1e+07 1000 10000 100000 le+06 1e+07 

degrees ot Ireedom degrees ol treedom 

(a) ground state energy (b) rjh^ i ly Q) 

Figure 5.4: The convergence curves of the ground state energy and ^(vE^, fi). 
Example 3: Fullerene Cqq. 

The ground state energy obtained by SIESTA is —341.340 a.u.. In our computations, we choose 
ti = [-30.0, 16.0] x [-23.0,22.0] x [-24.0,21.0] to be the computational domain. 

We can see the preservation of carbon- hydrogen bonds in Figure 15.51 which validates our calcu- 
lations. Figure 15.51 Figure 15.61 and Figure 15.71 show that more mesh points are placed around the 
atoms. 

The convergence curve of the ground state energy approximations is shown in the right of Figure 
15.81 from which we observe a convergence to —342.722 a.u., which is very close to the reference 
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Figure 5.5: Cqq: configuration and charge density on a sphere. 





Figure 5.6: Ceo'- charge density and mesh on a cross-section. 




Figure 5.7: Cqq: charge density and mesh on plane z = 0. 
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energy. The convergence curve of the a posteriori error estimator obtained by the quadratic finite 
element is shown in the left of Figure l5.8l from which we see that it reaches the optimal convergence 
rate. 



C60 C60 




100 1000 10000 100000 le+06 1e+07 °'\oo 1000 10000 100000 le+06 1e+07 

degrees of treedom degrees ol treedom 

(a) ground state energy (b) Jfe^h, Q) 

Figure 5.8: The convergence curves of the ground state energy and rjh(^h, 



6 Concluding remarks 

In this paper, we have studied the AFE approximations of Kohn-Sham equations. We have obtained 
the convergence and quasi-optimal complexity of the AFE approximations. We have also curried 
out some typical numerical simulations that support our theory. 

In our analysis of convergence rate and complexity of AFE approximations, for convenience, 
we have assumed that the numerical integration was exact and the nonlinear algebraic eigenvalue 
problem was exactly solved. Indeed, the same conclusion can be expected when the error resulting 
from the inexact solving of the nonlinear algebraic eigenvalue problem and the error coming from 
the inexact numerical integration are taken into account. 

Suppose that (A, $) € 0, the associated exact solution over mesh Th is (A^, $^), and the inexact 
numerical solution is (A/,, <!>/,). If the numerical errors resulting from the solution of (nonlinear) 
algebraic system and the numerical integration are small enough, say, satisfy 

\\$h - + |A h - A* | < r{ho)rft@ h , ft) 

with f(ho) <C 1 for ho <C 1, then we have from the following triangle inequality 

|A-A h | < lA-Afcl + IAfc-Ahl, 

and the similar perturbation arguments that the same convergence rate and quasi-optimal com- 
plexity can be derived. 

Finally, we emphasize that, in this paper, we have not given the convergence rate and complexity 
for the AFE approximations for the Lagrange multipliers A. Indeed, the related results for Lagrange 
multipliers are not so obvious, and we need do some more detail analysis, which will increase the 
length of this paper. We will consider it in our coming papers. 
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Appendix: A boundary value problem 

In this appendix, we shall provide some basic results for the AFE approximations of a model 
problem that was used in our previous analysis. Consider a homogeneous boundary value problem: 



£$ = F in Q, 
$ = on dtt, 



(A.l) 



where F = (fi)f =1 € (L 2 (Q)) N . Note that (|A3j) is equal to: Find $ € U such that 

a($,r) = (F,r) vre-H. (A.2) 

A standard finite element scheme for (|A.2[) is: Find € Vh satisfying 

a($, l; r) = (F,r) vrev ft . (A.3) 

Let T denote the class of all conforming refinements by bisections of 7o- For Tt € T and any 
T = {ji)iLi € VJu we define the element residual 1Z T (T) and the jump J e (r) by 

£r(T) = U + ^A^ mT€%, (A.4) 

/l 1 \ N 

Je(T) = ^2 V7i ' Tl ' ni + 2 V7i ' T2 ' "^J oneef/,, 

where e is the common face of elements t\ and T2 with unit outward normals n\ and , respectively. 
For r € 7ft, we define the local error indicator r)/j(r,T) by 

^(r,r) = ^||^ T (r)||^+ £ h e \\J e {T)\\l >e (A.5) 



and the oscillation osc/,(r,T) by 



6rc fe (r,r) =h T \\n T (T)-n T (T)\\ 0tT . (A.6) 
Given T' C Th, we define the error estimator fjh(T,T') and the oscillation osc/^r,? - ') by 
^(r,T') = J] rt(r,r) and 6Tc£(r,T') = ^ oTc^(r,r), 

rST' rGT' 

respectively. We see that a similar a posteriori error estimate to that for Poisson equation can be 
expected for (TjO]) (c.f. gTJ |H H3] ) . 

Theorem A.l. Let Q <E H be the solution of lA.fy) and <&h € Vh be the solution of l(A.3\) . Then 
there exist constants C\ , C2 and C3 > depending only on c a in (|4.1j) and 7* in ()2. 12|) such that 

||*-3fc||?^<<7i»7h(3fc,fi), (A.7) 
C 2 fjl($ h ,ty < \\$-$ h \\l n + C a 58cl($ h ,n). (A.8) 
An AFE algorithm for (IA.2I) is designed as follows (c.f. [T2"]): 
Algorithm A.l. 

1. Pick a given mesh To, and let k = 0. 
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2. Solve the system (j A.3|) on Tk to get discrete solution 

3. Compute local error indictors rjk(<&k,i~) for all r 6 7fc- 

4- Construct M.^ by Dorfler Strategy and parameter 6. 

5. Refine Tk to get a new conforming mesh Tk+i- 

6. Let k = k + 1 and go to 2. 

Using the similar arguments to those for scalar linear elliptic boundary value problem (see, e.g, 
[12]), we have the following result for Algorithm lA.il 

Theorem A. 2. If {$>k}keN * s a sequence of finite element solutions produced by Algorithm \A.1[ 
then there exist constants 7 > and £ £ (0, 1) depending only on the shape regularity 7* and the 
marking parameter 9, such that for any two consecutive iterations 

II* -**+i|li,n + T^+i < l a (ll*-4*lli,n + 7^(*fc.n))- 

Indeed, the constant 7 has the following form 

^(TT^W (A ' 9) 

with C* > depending on the regularity constant 7* and S £ (0, 1). 

For the distance between two nested solutions of (|A.3|) . we have (c.f. [T2] ) 

Theorem A. 3. Let <£># € Vh and $^ € Vh be solutions of IA.3\) respectively. IfTh is a refinement 
of Th by marked element A4h md refined elements TZ = TZ-j- H ^.j- h , then we have the following 
localized upper bound 



l*ff-*h|li,n<<5i 1"). 
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